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Abstract 

A new class of integro-partial differential equation models is derived for the prediction 
of granular flow dynamics. These models are obtained using a novel limiting averaging 
method (inspired by techniques employed in the derivation of infinite-dimensional 
dynamical systems models) on the Newtonian equations of motion of a many-particle 
system incorporating widely used inelastic particle-particle force formulas. By using 
Taylor series expansions, these models can be approximated by a system of partial 
differential equations of the Navier-Stokes type. The exact or approximate governing 
equations obtained are far from simple, but they are less complicated than most of 
the continuum models now being used to predict particle flow behavior. Solutions 
of the new models for granular flows down inclined planes and in vibrating beds 
are compared with known experimental and analytical results and good agreement is 
obtained. 

1 Introduction 

The last two decades have witnessed an intensification of research in granular flow dy- 
namics, in large measure spurred by a burgeoning array of engineering and industrial 
applications of particle technology. There are several features that make granular flow 
research attractive to engineers, mathematicians and scientists, among which are the fol- 
lowing: A need still exists to formulate the underlying principles of particle interactions 
in a completely satisfactory manner; there are as yet few if any definitive mathematical 
models that can reliably predict a wide range of granular flows; particle flow phenomena 
such as arching, surface waves and convection are still not entirely understood from a 
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mathematical or engineering perspective; there is a panoply of extremely complex non- 
linear dynamical behaviors exhibited in granular flow regimes that has not yet been fully 
analyzed and has severely tested or exceeded the capabilities of current experimental and 
computer technologies for accurate characterization; and the techniques and devices for 
optimizing certain features of particle flows are for the most part only understood on an ad 
hoc basis. In this paper we use some averaging and limiting ideas associated with infinite- 
dimensional dynamical systems theory to derive a new class of continuum mathematical 
models for granular flow that may be capable of predicting the dynamical characteris- 
tics of particle flows in a large variety of circumstances, and thereby help to make some 
progress in solving the many outstanding problems in this field. Our purpose is not to 
compete with the host of interesting models that try to incorporate as much of the physics 
of granular flows as possible, including (vibrational) energy equations. Rather, we aim to 
produce mathematical models that are relatively tractable and ignore just enough of the 
physics to still provide useful predictions of granular flow dynamics for a wide range of 
applications. 

Although there have been several partial successes in recent years, the state-of-the-art 
in mathematical modeling of granular flow phenomena pales in comparison to that of fluid 
mechanics where there is a universally accepted model - the Navier-Stokes equations - 
whose reliability has been tested and confirmed for over a century, and is considered in 
many quarters to be capable of apprehending even what may be the most elusive of all 
physical processes - fluid turbulence. Approaches based on continuum mechanics (trans- 
port theory) and kinetic theory (statistical mechanics) have been those most often used 
for obtaining mathematical models for particle flows in the form of systems of partial 
differential equations. Notable examples derived using these methods which have enjoyed 
some success in predicting granular flow dynamics may be found in An & Pierce Q, An- 
derson & Jackson Q, Farrel, Lun & Savage |^], Gardiner & Schaeffer |^, Goldshtein & 
Shapiro |10|, Jenike & Shield [12|, Jenkins & Savage |13], Jenkins & Richman |lj], Johnson, 
Nott & Jackson |l|l, LunJ^, Lun & Savage |l9|, Numan k Keller pi, Pasquarell pi. 
Pitman |24], Rajagopal |26], Richman [^^, Tsimring & Aranson [{j^ j. Savage ||2^, 3C], 
Savage & Jeffrey |^l|, Schaeffer |^2|, Schaeffer, Shearer & Pitman and Shen & Acker- 
mann |Q] (see also Fan & Zhu @ and Walton |3^] ) . Several of these models have proven to 
be rather effective in characterizing certain granular flows, for example in chutes and hop- 
pers with simple geometries, but they tend to be fairly complicated systems of nonlinear 
partial differential equations that are difficult to analyze and solve except by approxi- 
mate numerical methods, and the information they provide has barely made a dent in the 
host of practical problems associated with industrial uses of particle technology. There 
have also been a number of simple, idealized models formulated by neglecting a variety of 
physical factors, but these tend to miss many of the features of granular flow of interest in 
applications. Much work still remains in finding a really effective balance between math- 
ematical tractability and adherence to the underlying principles of physics in the models 
for a large class of granular flow phenomena. It is hoped that the models introduced here 
will provide a useful step in the direction of achieving such a balance. 

Granular flows have also been extensively studied using methods inspired by molecular 
dynamics research. The basic idea of the molecular dynamics approach is to use realistic 
models for interparticle forces, developed from both theory and empirical investigations, in 
a Newtonian dynamics context with a large number of particles (hundreds or thousands) 
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to determine the evolution of a particle flow configuration. Analytical means are of little 
use in solving the very high dimensional dynamical systems encountered in such an ap- 
proach, but some very sophisticated simulations, employing a variety of numerical solution 
techniques, have been devised for studying granular flows, such as those of Goldhirsch et 
al. Q, Lan & Rosato |16, 17|, McNamara & Luding |2^, Poeschel & Herrmann |25], Rosato 
et al. (2^, Swinney et al. |35] and Walton |^^. Alternative approaches based on cellular 
automata models and kinetic models of random walks in discrete lattices have also proven 
to be quite useful; see, for example, Baxter & Behringer ^ and Caram & Hong |^]. These 
and other simulations have proven to be so remarkably accurate in manifesting most of the 
complex aspects of particle flow behavior, that one is inescapably drawn to the conclusion 
that the formulation of a more concise and tractable mathematical representation of such 
simulations should greatly enhance our ability to analyze particle flow phenomena. 

It was this idea of finding more succinct ways of mathematically characterizing granular 
flow simulations for extremely large numbers of particles that served as the inspiration 
for the new models derived in this paper by computing limiting forms of the relevant 
Newtonian dynamical systems. To be more precise, we obtain systems of nonlinear partial 
differential equations - infinite-dimensional dynamical systems - for velocity fields of 
granular flows by using an averaging method together with the computation of a limit 
as the number of particles tends to infinity, followed by a Taylor series approximation. 
The approach employed is akin to the methods used to obtain limiting partial differential 
equations for systems of ordinary differential equations (as the size of systems tend toward 
infinity) in the theory of infinite-dimensional dynamical systems; for example, as when 
the Korteweg-de Vries equation is obtained as the "limit" of an infinite string of coupled 
nonlinear oscillators (cf. Tabor and Temam |37]). Our method leads to an infinite 
class of mathematical models of widely varying levels of complexity, depending on the 
form of the particle-particle force laws chosen and the order of the Taylor series expansions 
employed. Several of these models appear to enjoy certain advantages over existing models 
in terms of simplicity and ease of analysis, and they have the potential for providing a 
better developed mathematical understanding of granular flow phenomena. 

This paper is organized as follows: The particle-particle models, based on the Hertz- 
Mindlin theory and some empirical observation, that we shall employ for the granular 
flows under consideration are described in Section 2. Then, in Section 3 we develop 
the Newtonian differential equations of motion for the particle flow dynamics using the 
particle-particle force formulas introduced in Section 2, and we describe a decomposition 
of the forces into interparticle forces, body forces and transmitted forces. In Section 4 we 
delineate a limiting procedure on the Newtonian equations of motion of the granular flow, 
ignoring boundary contributions, that produces a system of integro-partial differential 
equations that models the velocity field of the particle flow. This class of mathematical 
models for particle flow dynamics is infinite and depends on the form and parameters of the 
particle-particle force formulas. We also show how our integro-differential equation models 
can yield an infinitude of partial differential equation approximations to the governing 
equations when the dynamical variables are approximated by Taylor series. A choice of 
interparticle parameters and order of Taylor series expansion that leads to a particularly 
simple system of partial differential equations for the particle velocities in a granular 
flow is also described in this section. We compute, in Section 5, an exact solution of 
the simple model of Section 6 for a fully developed flow through a vertical pipe with a 
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uniform circular cross-section, and in so doing give our first illustration of how to append 
relevant boundary conditions to our system of equations. The model of Section 4 and the 
introduction of appropriate auxiliary conditions necessary to model fully developed, two- 
dimensional granular flow down an inclined plane is treated in Section 6, and we compare 
our results with those obtained from other analytical and experimental studies of inclined 
plane flows. In Section 7 we develop the boundary conditions for flows in a vibrating bed 
and study numerically our model subject to these boundary conditions. We discuss our 
results and compare them with experimental studies. Finally, we conclude in Section 8 
with a discussion of the consequences of the work in this paper and possible directions for 
future research involving the new class of models. 



2 Interparticle forces 

In this section we shall describe the particle-particle force models that are the foundation 
upon which we construct our derivation of the governing equations of motion for the 
granular flows under consideration. We assume that the flow system is comprised of 
a large number N of identical inelastic spherical particles distributed throughout some 
region m at points x(*\ 1 < i < A^. The common radius of all the particles is a very 
small positive number that we denote by r, and the point x^*) corresponding to the ith 
particle is located at the center of the particle for all 1 < i < A^. 

We may select any particle, say the ith one, and suppose that the j'th particle at x^-'^ 
is near x^*). For convenience, we define 



r{ := x(j') - x(^) 

and to be the velocity of the jth particle relative to the ith particle; namely. 



(1) 



V X 



V X 



(2) 



Taking our cue from Hertz-Mindlin theory as supported by numerous experimental obser- 

] and |39|), we shall assume that 



12, 28 



vations and granular flow simulations (see 
the model for the force exerted on the ith particle by the jth particle is described as 
follows: is the sum of a (inelastic) normal force and a tangential force due to 
friction 



(3) 



where 



-X 



+ 1] 



V . r . 



(4) 



and 



(5) 



Here (•, •) is the standard inner product with induced norm || • || in and a, (3, 7 and 5 are 
positive exponents chosen according to the particular properties of the material particles; 
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among the most often used values are a = 1 or 3/2 (Hertzian), f3 = 1, 7 = 0, 1/3, 2/3 or 

3/2 and 5 = 1 or 2. The functions %, tp and rj are smooth (= C°°) on [0, 00) and have the 
following properties: 

x(T),^(T),r/(T) >0 for all r > 0; (6) 

X'(r),^'(r),r/'(r) < for all r>0 (^=d/dT)- (7) 

x(t) = V'(t) = 7?(r) = when r > 4r^; (8) 

x'(t) = V'(r) = r?'(r) = for < r < < 4r2; (9) 

r/(0) < x(0); (10) 

and 

r]{r) < x(t) for all r > 0. (11) 

Graphs of these functions are shown in Figure 1. As we arc going to ignore the rotational 
motion of the particles in our treatment, we shall assume that the tangential force is very 
small compared to the normal force and, more specifically, that ipi''') ^ vi^) all r such 
that iP{t) > 0. 

Figure 1: Force functions: a) y = x(x), h) y = r]{x), c) y = ip{x) 



The role of the function 77 is to represent an energy loss due to inelasticity in the 
restoring mode when the particles are separating after a collision. A caret over a vector u 
indicates the unit vector in the direction of u; i.e., u:=u/ ||u|| . The vector "d {^''i^ is the 
component of the relative velocity at the point of contact of a pair of particles obtained 
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by projecting onto the tangent plane of the ith particle at the point of contact. This 
vector can be written in the form 

^(v|) :=v|-(v^,f0f^ (12) 
The normal and tangential interparticle forces are depicted in Figure 2. 

Figure 2: Particle-particle forces: a) normal force, b) tangential force 



Summing over all particles in the granular flow system, we find that the total force 
exerted by all the particles on the ith particle is 

N N 

p-= E E (^+^)- (13) 

Note that our assumed particle-particle force models account for the geometry of the 
particles only with regard to the region where the force vanishes (its support) and the 
manner in which the tangential frictional component of force is defined. We observe 
that (13) can also be obtained from a specific force density field surrounding the zth 
particle with the force supplied by each grain equal to this specific density multiplied by 
the volume f^rr'^. We shall return to this point in the sequel when we compute limiting 
forms of the particle dynamical system. 

3 Newtonian equations of motion 

The motion of the particles in the granular flow field may be described by a system of 3A^ 
second-order, ordinary differential equations expressing Newton's second law of motion; 
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VIZ. 



x» = Fi :=Pi + Ti + Ei + Bi (1 < i < iV), 



(14) 



where ' = d/dt^ m is the mass of each of the identical particles, Pj is the force exerted 
on the ith particle by all particles in direct contact with it as described in the preceding 
section, Tj is the transmitted force on the ith particle exerted by connected arrays of 
particles in contact with one another that touch a particle in direct contact with the 
ith particle, Ej is the external or body force on the ith particle which is usually just the 
gravitational force (but may sometimes also include electromagnetic and other forces) and 
Bj is the boundary force exerted on the ith particle by fixed or motile boundaries in direct 
contact with it that delimit the region in space in which the particles can move. Observe 
that the variables on which each of the components of force depend can be described as 
follows: 



where the dependence on t in Bj occurs when the material boundary of the flow region 
moves with time such as in the case of particles moving in a vibrating container. 

We can write the Newtonian equations of motion in a more concise form by introducing 
the following vector notation: Define the vector X in R^^ to be 





is usually constant 



m 





Then (14) can be rewritten in vector form as 



X = *(X, X, t) := *p(X, X) + *t(X, X, t) + *e + *(>(X, X, t), 



(15) 



where 




is the interparticle force per unit mass, 



*T := m-HTi(X, X, t), . . . , TAr(X, X, t)) 



is the transmitted force per unit mass, 



*e := m-^ (Ei,...,Ejv) 



is the external force per unit mass and 



:= (Bi (x(i) , x(i) , i) , . . . , B^ (x(^) , x(^) , t)) 
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is the boundary force per unit mass. In theory, if initial values of X and X are specified, 
then (15) uniquely determines the ensuing motion of all the particles, at least for small 
values of |t| (see [^6| and |37]). However, for extremely large values of N the work required 
to integrate (15) - analytically, when in the rare cases that this is possible, or numerically 
otherwise - tends to be prohibitive. Thus it is desirable to find an infinite-dimensional 
limit in some sense for (15) as iV — > oo, presumably in the form of a partial differential 
equation, that may prove to be more amenable to analysis. This is precisely what we shall 
do in the next section. 



4 Limiting models 

We shall demonstrate how new models for granular flow phenomena can be obtained by 
applying a certain type of dynamical limit procedure to the Newtonian equations (15). 
The reader will no doubt notice at least a vague similarity between our method and the 
continuum limit used in the Fermi-Ulam-Pasta model to obtain the Korteweg-de Vries 
equation (cf. ||3^). To begin with, we restrict our attention to points in the interior of 
the granular flow region that are not directly affected by interaction with the boundary. 
Consequently, for the time being we ignore the boundary force contribution in (14) or (15); 
the boundary effects shall be considered in the sequel when we study specific boundary- 
value problems. 

Referring to (14), we assume that the body forces are exclusively gravitational and that 
the Cartesian coordinate system has been chosen so that the gravitational force acting on 
each particle has the form 

Ej = -mge, (16) 

where g is the acceleration of gravity and e is a unit vector in the opposite direction to 
the gravitational field. Now we select a point in the interior of the granular flow field 
corresponding to the ith particle (at time t) which is moving along a trajectory determined 
by the vector field 

x»=v(x«,t) (17) 

and the location of this particle at time t = 0. 

The interparticle forces on the ith particle at the point x^*) = x are given by (13). Since 
we are going to take a limit as the number of grains goes to infinity, we need to average or 
distribute these forces in a way that insures the existence of such a limit and is conducive 
to its computation. This can be done by smearing the particles into a continuum (assumed 
to be locally uniform) and considering the interparticle force field to be obtained from a 
specific force density field. To be more precise, we assume that each particle is surrounded 
by a specific force density field of the same form cP^=. Whence, the interparticle force on 
the ith particle can be written 



(18) 
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where AVj is the volume increment occupied by the jth particle and c > is a multiplica- 
tive factor with units volume~^ (associated with the geometry of the particles). This can 
be rewritten in the form 



(iV-l)-icco5;P. (x(^);v]) 



(19) 



where cq is a positive constant equal to the volume of the (compact) support of Pj and we 
have assumed that all the particles occupy volume increments of the same size. In (19) 

we plainly see the averaging aspect of this approach. Taking the limit as ^ oo in (19) 
[or equivalently as AVj- ^ in (18)] using standard results from integration theory, we 
obtain 



lim Pj 

AT— >oo 



c / P* dyidy2dy3 = c P* dy, 

JR3 JR3 



(20) 



where it follows from (3), (4), (5) and (12) that 



P* (y;v(x),v(x + y)) 

.l|2 



-X 



yf) (v(x + y)-v(x),y)||yf-i 



+i' (llyf ) lly| 



v(x + y) - v(x) - (v(x + y) - v(x), y) —2 

l|y|| 



y 

5-1 



(21) 



v(x + y) - v(x) - (v(x + y) - v(x), y) 



|y|| 



where y = iyi, y2,y3) represents the position vector measured from the reference point x 
that has been introduced to simplify the notation for . As for the transmitted force in the 
Newtonian equation (14), we make the standard assumption that in the continuum limit 
it can be represented by a gradient field, gradp, where the function p = p{x, t) is naturally 
called the pressure. The particle at x^*) is represented by a density field, p = p(x, t), with 
compact support. Hence the external force is 



pgdVj e, 

where Wi is a spherical (control) region centered at x^*) with (Lebesgue) measure AVi, 
and this converges to pge as N ^ 00 (4» AVi 0)- In the same spirit, the right-hand 
side of (14) is replaced by 



di 



I 

JWi 



pwdV, 



which upon applying the usual continuum limit converges to the density times the total 
(material) derivative of the velocity: 



dv 



dt 



k=l 



dxk 
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Upon combining all of the above computations, we obtain the following system of nonlinear 
integro-partial differential equations for the momentum balance of the particle flow in the 
interior of the region under consideration: 

d'v dv 1 f 

-FTT = ^ + Vk-K— = -9& + -gT^adp + K P*(y;v(x),v(x + y))d2/, (22) 
Dt dt dxk p Jr3 

where k := c/p and we have employed the Einstein summation convention. If p is constant 
and gradp is known a priori, then (22) together with appropriate initial and boundary 
data suffices to determine the velocity field. When the granular flow is compressible and 
gradp is known a priori, we have to add the continuity equation 

|^ + div(H = (23) 

to (22), and then by imposing additional auxiliary data the velocity field and density 
may be determined. Of course, in general, both p and p are unknown variables, in which 
case (22) and (23) are insufficient to determine p, p and v. One more equation must be 
added, and this may be accomplished by appending an energy equation to (22) and (23). 
The easiest way to do this is to obtain an equation of state of the granular flow medium 
that provides a relationship between the pressure and the density. If we make the same 
assumption as above (in particular, that the particles are uniformly and isotropically 
distributed locally), then by applying the same type of limit as iV — > co to the equations 
representing the kinetic energy of the Newtonian system (14), we obtain the equation of 
state of an ideal gas; namely 

p = Ap'^, (24) 

where A > and co > 1 are constants that are obtained from the properties of the granular 
flow medium. The same result can be derived by applying the standard thermodynamic 
limit of statistical mechanics in conjunction with the virial theorem. 

For certain purposes, including comparison with other continuum models for parti- 
cle flows, it is useful to replace (22) with an approximate partial differential equation. 
Although, it should be pointed out that, mathematically speaking, (22) enjoys certain 
inherent advantages over such partial differential equation models. In particular, as we 
shall demonstrate in a forthcoming paper, solutions of the system with (22) exhibit con- 
siderably more regularity than the pure differential equation models that we shall discuss 
in the sequel. 

In order to approximate (22) by a system of mth order partial differential equations, 
we may use the following Taylor series expansion of order m : 

v(x + y)-v(x)^^--J(x)y^ (25) 

k=l 

The positive integer m is at our disposal, and it is plausible to assume that the larger we 
choose m, the more accurate the approximation. Substituting (25) in (22), we obtain the 
system of mth order, nonlinear partial differential cqiiations as an approximate model for 
the momentum balance of the granular flow field given by 

dv Oy ~ ^ 1 f dv d™-v\ 
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where T, a function that does not depend expUcitly on x, is defined by 




Hence we have infinitely many possible partial differential equation models for granular 
flow corresponding to the choices of the functions x, rj and of parameters a, (3, 7 and 5, 
and the order m of the Taylor series approximation. This leads to a very natural question: 
What order of Taylor series approximation in (26) should be used for a given application? 
As we shall show in the sequel, m = 2 works rather well for tube, inclined plane and 
vibrating bed flows. However, it will probably be necessary to consider several choices in 
other applications and determine an acceptable order of approximation on a case-by-case 
basis, where an educated guess is made based upon known properties of the flow. 



A simple flow model 

Depending on the choice of parameters and the order, the equation (26) can range from 
relatively simple to quite complicated. In this section we make a choice of parameters and 
order that leads to a rather simple yet ostensibly realistic model for the velocity field of 
a granular flow. Specifically, we choose q; = /3 = 1, 7 = 0, 5 = 1 and m = 2. Then (26) 
takes the form 



dt dxk 



- 1 

-ge -\ — grad 

P 



p + K 



-X 



^ 1 



fe=i 



which upon integration using spherical coordinates simplifies to 

d\ 9v ^1 , ^ ^ , / T N 

-- — h Vk — = "I — grad p + Av + A grad ( div v) , 

ot OXk P 



where 



u := 



2'KK. 

T5~ 



/ H 

Jo 



2\ 2 

S )s 



4V(s^)] sUs 



and 



A := 



Attk 
T5~ 



[r/(s2)s2 - V(s^)] sUs 



(27) 



depend only on the density and the particle-particle force functions described in Section 2 
and may be assumed to be constants in many applications. 
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Observe that (27) is essentially just the momentum part of the Navier-Stokes equations 



(cf. [p3| , 36] and ^7^). There are at least two interesting inferences that may be drawn 



from this result: Firstly, it provides a partial confirmation of the validity of our integro- 
partial differential equation model as a predictive tool for granular flow. Secondly, it lends 
support to the contention that the Navier-Stokes equations are a good model for granular 
flow behavior obtained from simulation of the Newtonian equations of motion. 



5 Flow through a tube: an exact solution 

In this section we obtain an exact solution of the approximate model (27), (23), (24) 
subject to appropriate boundary conditions, for the case of fully developed (steady-state) 
granular flow, under the action of gravity, through a vertical circular cylindrical pipe 
illustrated in Figure 3. We assume that the density and pressure are constant, hence it 
suffices to solve (26) subject to some boundary conditions. 



Figure 3: Flow through a tube 



Under the circumstances, it is convenient to recast (27) in terms of standard cylindrical 
coordinates (r, 6, z) with corresponding velocity components {u, v, w), where u is the radial, 
V the azimuthal and w is the vertical(axial) component of the flow velocity. The system 
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assumes the following form with respect to cylindrical coordinates: 



1 d'^u 



r dr \ dr J dO"^ 



+ 



dv dv V dv dv uv 
ot or r o(f oz r 



Id ( dv\ 
-T^ r— + 



r dr \ dr 



1 d'^v 



d'^v 
d^ 



+ A— 
dr 



Xd_ 
^ rde 



1 d{ru) 
r dr 



1 dv dw 
r do dz 



1 d{ru) 
r dr 



^ 1 9f ^ dw 
r dd dz 



(28) 



dw dw V dw dw 
dt dr r dd dz 



-g + v 



ld_ 

r dr 



dw\ 
dr J r 



1 d^ 



w 



+ 



92 



w 



2 QQ2 



dz 



1 d{ru) 
r dr 



w 



1 dv d 

r do dz 



where g is the acceleration of gravity. We assume that the pipe has radius R> Q and that 
its length is so great that the domain of the granular flow can be represented in idealized 
form as 

n := {{r,0,z) ■.0<r<R}. 

Now we deal with the task of appending appropriate auxiliary data to (28) on ft. As we 
are seeking a steady-state solution, we assume that the velocity is independent of the time t. 
There remains the question of realistic auxiliary data on the boundary dQ. Of course, u <0 
on dfi is required by the geometry of the pipe (assuming it is rigid and impenetrable). 
Over time, one may reasonably expect the radial and azimuthal fluctuations in velocity 
along the inside surface of the pipe to cease, so we shall assume that both u and v vanish 
on dQ. As for the axial velocity along dU: the motion of a particle in contact with dft is 
that of free fall with a resisting force due to friction. This suggests that there is a constant 
limiting (or terminal) velocity along the wall of the pipe (that is achieved in the long-term 
flow configuration), so it is reasonable to assume that w is a negative constant along dfl. 
In summary, we take the auxiliary data for (28) in f2 to be 

du dv dw . ^ 

m il, 

(29) 



dt dt dt 



u = v = 



and 



W = —Wn 



on 



dn, 



where Woo is a positive constant. 

In view of the boundary conditions, it makes sense to seek a solution of (28) (29) with 
u = V = and w = (p{r). Then the first two equations of (27) are trivially satisfied and 
the third equation yields 



v d 
r dr 



dip 
dr 



9- 



(30) 



Integrating (30), we obtain 

9 1 1 
= — r -I- ci logr + C2, 
4i/ 
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where ci and C2 are constants of integration. The solution should be regular at r = 0, so 
we must set ci = 0. Then (29) leads to the following solution: 



It is easy to check that (31) satisfies (28)-(29). 

6 Inclined plane flow model 

The fully developed flow of particles down a two-dimensional inclined plane will be studied 
in this section using the governing equations (27), (23), (24). As inclined plane flow has 

been extensively investigated (sec, for example [2] and [10]), we shall have an opportunity 
to compare the predictions based upon the simple approximate model with the results 
obtained by other researchers, thereby further testing the effectiveness of our approach. 
Figure 4 depicts the flow geometry for a plane inclined at an angle of to the horizontal. 

Figure 4: Inclined plane flow 



It is convenient to use a Cartesian coordinate system with x measured down along the 
surface of the inclined plane and y-axis normal to the plane and pointing into the flowing 
layer of granular material. Here u represents the component of the flow velocity along the 
X-axis and h the depth of the flowing layer. We assume that the density is constant and 
that the pressure gradient exactly balances the gravitational force normal to the inclined 
plane throughout the flowing layer. Therefore, it suffices to solve (27) along with the 



u = V = Q 



and 




(31) 
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necessary boundary conditions. Of course, the figure embodies the usual assumption that 
the granular flow is essentially two-dimensional. 

A clockwise rotation of 9 of the coordinate system and a balancing of the gravitational 
and reaction forces normal to the inclined plane yields the following pair of equations for 
the granular flow: 

du du du ■ Q f d'^u\ ^ / d'^u d'^v \ 

dt ^ dx ^ dy ^ ™ ^ \ dx'^ dy"^ ) \ dx^ dxdy / 

dv dv dv / d'^v d'^v\ ^/ d'^u d'^v\ 
dt dx dy \dx'^ dy'^ J \dxdy dy"^ J ' 

where V is the y-component of the velocity of the granular flow. Since the flow is taken 
to be fully developed (steady-state), it is reasonable to assume that both u and v are 
independent of the time t. It is also sensible to presuppose that the y-component of the 
velocity vanishes identically and that u is a function of y only. With these assumptions (33) 
reduces to the trivial equation = and we are left only with the simple ordinary 
differential equation 



(fu _ g 
dy'^ V 



2 - sin 6* (34) 



representing (32). 

Appropriate auxiliary data for (34) are the free-boundary condition along the free 
surface representing the interface between the flowing particles and the air that defines 
the depth of the flowing layer h as the smallest number satisfying 

fin 

-{h) = {h> 0), (35) 

and a slip condition along the inclined plane granular material interface 

^(0) = kuoo, (36) 

where Uqq is cl limiting velocity along the surface of the plane. The component u^o may 
be the result of a partial balance between the frictional properties of the plane and the 
particles and the gravitational component of the force in the x-dircction or a combination 
of gravitational and frictional effects and some constant mass flow rate supplied to the sys- 
tem. Here k is some nonnegative constant connected with the nature of the shearing stress 
in the flowing layer adjacent to the plane that is related to the frictional characteristics 
of the plane and particles and the dynamical state of the system. 
Integrating (34) twice using (36), we obtain the solution 

u = u{y) = -(^^]y'' + u^{ky+l). (37) 



Whence we determine the depth of the flowing layer by substitution of (35) in (37) ; namely, 

h = — — , (38) 
gsm9 
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for ^ > 0. A typical velocity profile is shown in Figure 5. 
Figure 5: Velocity profiles for inclined plane flow 
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The extremely simple nature of the solution (37) obtained from the governing equa- 
tion (27) notwithstanding, it compares rather well with observations from experimental 
studies and the predictions from more complicated flow models (cf. |^ and [^]). For 
example, the form of the velocity profiles illustrated in Fig. 5 is qualitatively similar to 
those measured in experiments and derived from more comprehensive constitutive equa- 
tions. Moreover, unlike some fairly popular models, (38) shows that our approach predicts 
a decrease in the depth of the flowing layer with increasing inclination angle of the plane - 
a property that is consistent with experimental observations. 

7 Vibrating bed model 

In this section we shall apply our model (27), (23), (24) to the study of granular flows in 
a two-dimensional vibrating bed. Namely, we consider the motion of a very large number 
of particles in a rectangular container in the plane with fixed vertical side walls and a 
horizontal bottom that is oscillating periodically in the vertical (X2) direction. The only 
body force is a gravitational force in the negative {X2) direction and the interstitial and 
surrounding medium is air which we assume has no effect on the granular flow. 
At i = the particles are contained in the following region: 

Kq := {{xi,X2) G R2 : |xi| < a, X2 > 0}, (39) 

where a > is half of the width of the container. Then, the container is subject to a 
vertical oscillation of the form asin(a;t) that is illustrated in Figure 6. 
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Figure 6: Granular material in a vibrating container 



We shall use boundary conditions at the walls similar to those employed in the previous 
section. Namely, we assume that the normal component of the particle velocity near the 
wall is equal to the normal component of the wall velocity. As for the tangential direction, 
we use the equation (36) in the following form 



^ = -kvT. (40) 



where vt denotes the relative tangential component of velocity between the particle and 

the wall, d/dn is the partial derivative in the outer normal direction. The equation (40) 
as well as (36) represents a type of balance law between the interparticle and particle-wall 
friction forces which has been used by other reseachers. The (constant) coefficient k > 
is a measure of the boundary friction that we shall call the wall friction coefficient. It 
would be most natural to use the particle size as the characteristic length for the non- 
dimensionalization of k, but it tends to zero in the continuum model limit. Therefore we 
use for this purpose the space step of the numerical integration scheme. 

In summary then, we take the following as the governing equations plus the initial and 
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boundary conditions for the granular flow in the planar vibrating bed: 

dvi dvi dvi ( d'^vi d'^vi \ , / d'^vi d'^vo 
— - + ?;i + vo = u ^ A ^ + A i- A 



dt ^ dx\ ^ dx2 \ dx\ ~'~ dx^ ) ^ \ dx\ ^ dxidx2 
dv2 dv2 dv2 f d'^V2 d'^V2\ , / d'^vi d'^V2\ 

1- 7;i 1- l)n = u ^ -I ^ -I- A 1 ^ 



dt dxi ^"^8x2 ^ V (^'^\ ^ 9x2 ) ^ ^ \ dxidx2 ~^ dx^ J 



(41) 



in E := {{x,t) : x £ ^It, t > 0}; 

vi(xi,a:2,0) = t;i(xi,a;2), ^^2(2:1, X2, 0) = ^^(xi, X2) at t = (42) 

for all a; G = {x G : < a, < X2 < h}, where the functions Vi{xi,X2) and 
V2{xi,X2) determine the initial velocity distribution; 

— — = —kvi, V2 = OLOcosiut) (43) 
0x2 

for all particles on the bottom, {{xi,X2) \xi \ < a, X2 = a sin(a;t)}, of the bed when t > 0; 

VI =0, ^ = -kv2 (44) 
0x1 

for all particles on the left wall, {{xi,X2) ■ xi = —a, X2 > asin(a;t)}, of the container 
when t > 0; 

^'i = 0, ^ = -kv2 (45) 
0x1 

for all particles on the right wall, {{xi,X2) xi = a, X2 > asin(a;t)}, of the container 
when t > 0; and 

I- 

for all particles on the free-boundary, consisting of all points in dQ,t\dKt, when i > 0. 

Now we shall consider some numerical solutions to the model (41) subject to the initial 
and the boundary conditions (42)-(46). To simplify our analysis, we shall ignore the 
effects of surface waves and free-boundary components at the bottom of the container. 

In order to avoid difficulties with vibrating bed boundary conditions like V2{x,a sin(a;i)) 
= auj cos{ujt) at the bottom of the bed, we shall write our equations in the vibrating system 
of coordinates: 

x = x*, y = y* + asm{Lot), t = t*. (47) 

The operators of partial differentiation in this frame are 

d _ d 9 _ 9 d _ d ( *\ ^ /'/1Q^ 

dx dx* ' dy dy* ' dt dt* dy* 

We have to include also into the system (41) the inertia! force term proportional to 
au'^ sin(a;t) and directed along the y-axis. The governing system of equations in the 
"starred" system takes the following form (we have dropped the index '*') 

Vl,t = t^ivi^xx + Vi^yy) + aUJ COs{ujt)vi^y " a{viVi^x + V2Vl^y) + X{vi^xx + V2,xy), 
V2,t = aLO^sin{ujt) + y{v2,xx + V2,yy) + aUJ COs{ujt)v2,y - Ci{viV2,x + V2V2,y) (49) 
+Kvi,xy + V2,yy), 
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where 0<x<2, 0<y<2, and the boundary conditions can be written as 
vi{0,y,t)=0, vi{2,y,t)=0, V2{x,0,t)=0, V2{x,2,t)=0, 
-^(x,0,t) + /c?;i(x,0,t) =0, ^(x,2,t) = 0, ^^^^ 

^(0, y, t) + kv2{0, y, t) = 0, ^(2, y, t) + kv2{2, y, t) = 0. 

We use an exphcit finite difference scheme for solving the system (49)-(50) with the 
spatial step Ax = Ay = 0.001 and the time step At = 10~^. Estimates show that such a 
small time step is needed to satisfy the stability condition for the explicit finite-difference 
scheme. 

We investigated the system (49)-(50) using both multi- vortex and random initial condi- 
tions. For certain ranges of the parameters (corresponding to the particle-particle forces) 
the motion of the system starting with a multi- vortex configuration changes rapidly into a 
pair of vortices that persists for a long time (relative to the period of the forced oscillations). 
The centers of this "stable" vortex pair oscillate with small amplitude synchronistically 
with the forced oscillations. In some cases this vortex pair evolved into a single vortex over 
a very large time period. Increase of the constant results in a corresponding increase of 
the particle-particle friction and leads to damping of the vorticity (see Figure 7) . 

When the motion starts from a random initial velocity distribution we also observed 
the "stable" vortex type of motion and bifurcation between different types of relatively 
stable patterns (Figure 8). Some of the values used for the control parameters were a; = 2, 
a = 1, V = 0.3, A = 1.0 in the case of the vortex type of motion, and a; = 3, a = 1, 
v = 1.0, A = 0.5 in the case of a mixing motion. These types of particle dynamics are in 
agreement with experimental observations and computer simulation results |1T7| , We 
note also that Hayakawa and Hong |11] obtained similar results from numerical solutions 



of their models but they assumed no-slip boundary conditions at the walls which are not 
physically realistic for vibrating bed granular flows. 

Similar types of the flow behavior were obtained by Bourzutschky and Miller Q for 
their Navier-Stokes models. By using negative slip boundary conditions in numerical 
experiments corresponding to granular flows with a high mobility boundary layer, they 
obtained experimentally observable vortex type solutions. Unlike us, they did not obtain 
convective flow behavior coinciding with experimentally observed results for possible val- 
ues of the wall friction coefficient. A possible explanation for this discrepancy between 
their findings and ours may be the fact that we included the gravity force directly in our 
model and they did not. 

As mentioned above, we have suppressed the free-boundary conditions that occur in 
an actual vibrating bed in our numerical experiments. This has been done to simplify the 
numerical solution of the problem, since incorporation of the free-boundaries significantly 
complicates the problem and for us is still in the developmental stage. Preliminary results 
indicate that the addition of free-boundaries will still result in the appearance of "stable" 
convective vortices, and we plan to demonstrate this in a forthcoming paper. Apparently, 
the frictional effects of the walls is the primary mechanism in the generation of convective 
rolls. For now then, our analysis of the vortices must be considered to be of a local rather 
than a global nature. 
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Figure 7: Motion of particles starting with four-vortex initial configuration (numerical 
solutions) 
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Figure 8: Motion of particles starting with random initial velocity distribution (nu- 
merical solutions) 
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8 Concluding Remarks 

Starting with well-estabhshed representations for particle-particle normal and tangential 
frictional forces (based on sound theoretical principles and a large body of experimental 
observations), we derived a new class of integro-partial differential equations to describe 
the velocity field in the granular flow of rough, inelastic particles. These granular flow 
models were obtained by taking a dynamical limit as TV — > oo of the Newtonian system 
of differential equations of motion of an A^-particle array using integral averages of an 
assumed uniform distribution of particles comprising the flow field. Then by employing 
Taylor series expansions of key variables of the flow field, we were able to obtain an 
infinite collection of approximations of the model equations in the form of a system of three 
nonlinear partial differential equations for the velocity components of the granular flow. 
The simplest of these approximate models, obtained by retaining only the first two terms 
in the Taylor expansions, is a system of equations that is significantly less complicated than 
most of the continuum models currently being used to investigate particle flow dynamics. 

Our models, and especially the simplest of the approximations, certainly do not incor- 
porate as much of the physics involved in granular flows as do the more comprehensive 
partial differential equation models, yet they appear to be quite promising instrumentali- 
ties for the prediction of particle flow behavior. A good indication of this is the results of 
our application of the simplest model to granular flow through a vertical tube, fully devel- 
oped flow down an inclined plane and flow in a vibrating bed which produced relatively 
simple solutions that compared remarkably well, in a qualitative sense, with experimental 
observations and the predictions of more complete models. This suggests that, in spite of 
the simplifying assumptions we used in the derivation, these models may be capable of ac- 
curately predicting dynamical properties of a wider range of granular flow configurations 
than one might imagine. And that they certainly warrant further investigation and test- 
ing. Moreover, the new models are far more amenable to analysis (particularly from the 
viewpoint of infinite-dimensional dynamical systems theory) than the majority of govern- 
ing equations in the literature. Therefore it is quite possible that one may be able to 
apprehend important new insights into several elusive granular flow phenomena from a 
more penetrating mathematical investigation of the properties of the new equations. 

In the near future we plan to undertake an intensive analytical and computational 
study of the models introduced in this paper. For example, we shall use dynamical systems 
theory to identify and analyze such phenomena as inertial manifolds, bifurcations, strange 
attractors and regimes of spatio-temporal chaos that will then be correlated to a variety of 
complex granular flow behaviors. In addition, it should be useful to develop and implement 
algorithms for the approximate numerical solution of the models, and then compare the 
results obtained with those from simulations,experimental studies and other governing 
equations. We shall begin this program of investigation by conducting a more thorough 
analysis of vibrating bed flows and also studying granular flow in hoppers. 
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